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Abstract 

We compare the predictions of linearized theory for the radiation produced 
in the collapse of a spherically symmetric scalar field with a full numerical in- 
tegration of the Einstein equations. We find power-law tails and quasinormal 
ringing remarkably similar to predictions of linearized theory even in cases 
where nonlinearities are crucial. We also show that power-law tails develop 
even when the collapsing scalar field fails to produce a black hole. 
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I. INTRODUCTION 



Linearized perturbation theory has been the main analytical - and until comparatively 
recently, numerical - tool for analyzing nonspherical gravitational collapse. The complexity 
of the problem has usually made this approach necessary, and it has been assumed until 
recently that the approach was sufficient. Recently, however, Gomez and Winicour |I| have 
focused attention on the extent to which these results are even qualitatively representative 
of the late stages of collapse. 

In the picture given by linear perturbation theory of the late stages of collapse, there 
are two features which are noteworthy. One is the development of "quasinormal (QN) 
oscillations," damped oscillations at complex frequencies characteristic of the mass of the 
black hole background. The second feature is the decay in time t of perturbations as 
the "power-law tails." There are good reasons to examine more carefully whether these 
features also appear in the fully nonlinear case. The arguments for the QN oscillations and 
for the tails are somewhat different, and should be considered separately. 

According to linearized theory the QN frequencies are fixed complex numbers multiplied 
by the inverse of the mass of the black hole background. (We use here and throughout 
units in which c = G = 1.) It seems reasonable that the phenomenon of QN ringing will 
be a feature of nonlinear collapse. One argument is that some numerical investigations of 
solutions of the fully nonlinear equations have shown QN ringing to be common [||. Secondly, 
the idea of QN ringing seems "robust." It is a natural frequency associated with a radiative 
boundary condition, and can occur in many different radiative systems. If QN ringing is 
found, to what black hole mass does it correspond? The QN oscillations themselves carry 
energy and may change the meaning of the mass. A reasonable guess, at least, can be made 
that the QN frequency evolves somewhat during the collapse. 

The situation for the power-law tails is quite different. These tails are not familiar or 
common phenomena. The explanation of their existence can be given in two very different 
ways: (i) They can be viewed as the result of the scattering of gravitational waves off the 
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"effective curvature potential" of the black hole spacetime 0, or (ii) they can be associated 
with the branch cut in the Green function for the wave propagation problem . Both argu- 
ments leave open the possibility that the tails are idiosyncrasies of the linear approximation. 

If QN oscillations or tails are missing from a fully nonlinear collapse, or if there is any 
significant new qualitative feature, the result might be to undermine confidence in the picture 
of collapse given to us by the analysis of linear perturbations of black hole backgrounds. 
Gomez and Winicour DH have addressed this question with numerical studies of the collapse 
of a spherically symmetric scalar field due to its own gravitational pull. Since the spherically 
symmetric problem involves only a 1+1 hyperbolic system, it is enormously easier to solve 
numerically than the problem of nonspherical collapse. 

What is more, the problem is a wonderful testing ground for comparing nonlinear results 
and the predictions of linearized theory. In linearized perturbation theory the evolution of 
a scalar field is governed by essentially the same mathematics that governs the dynamics 
of nonspherical perturbations. In particular, perturbation theory makes very specific pre- 
dictions about QN ringing and tails for perturbative spherically symmetric scalar fields. It 
is therefore of great interest that in their initial numerical studies of scalar field collapse, 
Gomez and Winicour have seen neither QN ringing nor tails. 

In this paper we will study the fully nonlinear evolution of a scalar field minimally coupled 
to general relativity. We consider first the evolution of a spherically collapsing scalar field; 
in addition, we consider nonspherical perturbations of this spacetime. We establish that 
the QN frequencies and the power-law tails of the numerical simulations are in remarkable 
agreement with the predictions of linearized theory when a black hole develops. If a black 
hole does not develop, and all energy eventually radiates away to infinity, we find that power- 
law tails still form. The existence of tails, but not QN oscillations, when holes do not form 
agrees with the analysis presented in the companion paper, hereafter referred to as Paper I. 

The organization of this paper is as follows. In section II we describe the coordinate sys- 
tem and the version of the field equations we use. In section III we describe our discretization 
and discuss the numerical error. In section IV we study the collapse of a scalar field for var- 
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ious initial configurations. In section V we study the evolution of multipole moments of test 
fields on the collapsing background. In all cases comparisons with the linearized results are 
made. We end in section VI with a brief summary and with conclusions. 



II. FIELD EQUATIONS AND ALGORITHM 



We study a spherically symmetric scalar field (p satisfying the minimally coupled equation 



^'1, = 0- 



(1) 



To describe the spherically symmetric spacetime we use the line element 



(is^ = —g{u, ^)g{u., r) du^ — 2g{u, r) du dr — dVt^, 



(2) 



in which u is a retarded time null coordinate. Regularity at the center requires that g = g 
at r = 0. The remaining coordinate freedom is fixed by the choice that u be the proper time 
on the r = central world line, ot g = g = 1 a.t r = 0. We introduce the auxiliary field ip by 



1 r 
= - ip dr. 
r Jo 



In terms of this variable, the wave equation (0) for takes the form 
The metric coefficients g and g are determined from (j) by 



g = exp 



4:71 I — r dr 
\ or 



(3) 



(4) 



(5) 



g = - f gdr. 
r Jo 



(6) 



Because of the spherical symmetry there is no independent gravitational degree of freedom. 
The operator D in differentiates along ingoing null characteristics, and information 
is propagated along outgoing null characteristics (the u = constant lines) by @. The 



integration scheme of equations and can therefore be related to a fully characteristic 
coordinate system {u = "retarded time", v = "advanced time,") in which the metric takes 
the form 

ds^ = —f{u^v) dudv + r'^{u^v) dVL^^ (7) 

and in which D oc d/dv\u- Here the coordinate v is fixed only up to arbitrary transformations 
V — >■ v{v). Since the coordinate is only a label on the ingoing null geodesies it can be assigned 
values f = 1, 2, 3... on our numerical grid. 

Our algorithm works on a characteristic grid made up of lines of constant u and f , with 
the radial "coordinate" r treated as a metric function r{u,v) which is evolved as 

Dr = -g/2. (8) 

As initial data for the algorithm we take ip{v) and r{v) on an initial outgoing null cone u = uq. 
These are equivalent to the choice of null data 0(r) and the (arbitrary) numerical choice of 
the initial position in r of each ingoing null lines of the grid (i.e., each line of constant v). 
The algorithm, which closely resembles that of Goldwirth and Piran proceeds as follows: 
We start by using (^ to obtain and by using (|^) and (||) to obtain g and g as functions 
of V along u = Uq. (The integrations over r are discretized as summations over v.) We next 
choose a value of the "time step" Am, and we use (H) to obtain ip, and (||) to obtain r, at 
grid points along u = Uq + Au. The process is then repeated starting with the new outgoing 
characteristic u = Uo + Au. On each such cycle of integrations the step size Am is chosen so 
that 

\r{v, u + Au) — r(t>, u)\ < \r{v + l,u) — r(t>, u)\ (9) 

for all f = 1, 2, 3, ... This has been found to be a useful rule of thumb. Since our integration 
scheme proceeds along characteristics, there can be no violation of causality, and hence there 
is no formal "Courant condition" to satisfy. 

Figure 1 shows our uv grid embedded in the well-known conformal diagram of spherical 
collapse. Null lines u = const, and v = const, are at 45 degrees. Infinity has been brought 
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to a finite distance and each point, apart from the hne r — 0, corresponds to a 2-sphere of 
surface Airr'^. This diagram shows where our null data are set and how far they are evolved. 
The upper left and right sides of our coordinate patch, though at finite distances, can be 
taken as approximations of the horizon 7i+ and future null infinity scri+. Lines of constant 
r go from past to future timelike infinity Those with r < 2M/, where Mf is the final mass 
of the black hole, cross the event horizon, while r = 2Mf approaches it asymptotically This 
gives one numerical method for finding the final black hole mass. Another is to take the 
limit of the Bondi mass Mb{u) at late retarded times w — > oo, where 



Yet another spacetime coordinate will be useful in comparing our results to analytic 
predictions, although it plays no active role in our algorithm. This additional coordinate is 
defined as the proper time along an r = const, trajectory, or 



We shall also use Bondi time, tB = t{u, oo), which is the retarded time coordinate that 
agrees with time at infinity for constant r. In an asymptotically flat spacetime the large r 
limit of t{r) is given by dts = lim^^oo g{''^, u) du. 

Our code has limited integration time precisely because it is characteristic. By definition 
u is the proper time of an observer at r = 0. The event horizon starts spreading out from 
r = at a finite value of u, say Uh. In a collapse g{oo,u) becomes infinite as u ^ Uh. In 
our algorithm this means that the stepsize Au decreases rapidly, while Ats remains about 
constant. The numerical approach to the horizon is stopped eventually by an overflow of g or 
underflow of Au. The situation is best illustrated in Fig. 2. The u = const, lines of our grid 
are squeezed together against the horizon. The problem of over/underflow, and of numerical 
instability (which usually develops earlier) is shared by all codes which avoid singularities 
by staying outside apparent horizons. A possible cure would be a slicing which does cross 
the horizon, combined with simply discarding parts of each slice inside the horizon. If the 




(10) 




(11) 
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surface beyond which time shces are discarded eats into the remaining part of the grid with 
the speed of hght or faster, no boundary condition is required on it [|^. Such a procedure 
is not apphcable to a characteristic grid, where the "time" shce is aheady an outgoing null 
cone and cannot be intersected by one. 

III. DISCRETIZATION AND ERROR ANALYSIS 

As we have just seen, our grid is highly nonuniform in m if a horizon forms. In that case 
it becomes also highly nonuniform in r. Ingoing null geodesies pile up at r = 2Mj, even if 
V is chosen uniform in r on the original slice. Truncation error is reduced, of course, if the 
grid spacing is made finer. We have tried to produce a code that is second-order accurate 
under a uniform grid rescaling. We denote here the relative size of the grid spacing by /i; 
a reduction by a factor 2 oi h means that the spacing in r is reduced by 2 and, due to the 
Courant-like condition (^), the spacing Au is also reduced by 2. 

On the top level we treat each of equations (^) and (^) as an ordinary differential equation 
in u at each fixed v. These are solved by the second-order Runge-Kutta method. The 
calculation of the "coefficients" 0, r, g and g in these "ODEs" is nontrivial however, and 
couples the equations at different values of v. They are given by the definite integrals (H), 
(I) and (H), which must be evaluated in that order. The integrals are discretized by the 
trapezoidal rule, which is again formally second-order accurate. 

In principle our code should be second-order. That is, for the finite difference approach 
we use, the error at a given spacetime point should vary as /i^. We did not find this second- 
order convergence one would expect naively from the code, but did find convergence better 
than first order. One reason for the failure of a simple error analysis is the ambiguity about 
the measure of error. As the dynamical range of cf) and the other fields is very large in a 
collapse spacetime, there is little point in taking an /2, or ^oo norm of the error in the fields 
over a set of spacetime points covering all of the integration region. The regions where the 
fields are large would dominate the integrated error. What is relevant to our confidence in 
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the results is the relative error, in particular in the regions where we are looking for tails 
and where 4> is small. 

Because of this we have looked at small spacetime regions, taking the I2 norm over only 
a few neighboring points. For all these regions we found power-law dependence of errors to 
apply over a change in grid size by a factor of 16. Roughly speaking, error scaled as h'^ , with 

= 1 for small r or small ts, and increasing up to 1.6 for both r and ts large. By small 
r we mean the very narrow shell in r, just outside r = 2M/, in which remains large at 
all times and where its gradient increases without bound. By small we mean the region 
where the bulk of outgoing radiation passes. 

Figures 3 and 4 display the error in of 0(r, i^), for sections of the hues r — 10 and 
ts — 50. In these figures the initial data is a Gaussian 4){r,tB — 0) with center 1.0, width 
0.1 and amplitude 0.06. (See the next section for a general discussion of initial data.) As 
there are no useful analytical solutions available, we had to use a numerical solution as the 
reference solution. We chose a solution generated with an initial grid spacing of Ar = 1/320, 
and compared to it grid spacings of 1/20, 1/40, 1/60, 1/80, 1/120, 1/160 and 1/200. (In all 
of these, however, the grid spacing was smoothly reduced by another factor of 16 for grid 
values of < r < 1, as we shall explain below.) The values of for the lower- resolution runs 
were interpolated linearly to the values of r and of the grid points of the highest-resolution 
run. The difference of (p between a run and the reference run was squared, summed and the 
square root taken. 

Figures 5 and 6 establish visually that the code converges, by showing sections of 0(r = 
10, ts) for different grid spacing. Figure 5 shows the region dominated by QN ringing, and 
Fig. 6 a region dominated by the power-law tail. 

It is perhaps surprising that the error varies as even if N varies over spacetime and 
does not have the value 2 that one might suppose. The fact, however, that this power is 
closer to 1 both near the origin and in regions of large 0, and closer to 2 elsewhere, seems 
to be an indication that the error is largely due to our handling of the r = boundary 
conditions. 
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We were alerted by reference to the possibility that the main source of error may 
be the boundary conditions g = g and (p = a.t r = 0. It can be seen in Fig. 1 that one 
after another, the v = const, grid lines cross r = into unphysical negative values of r, 
at which point they are dropped from the algorithm. We have numerically implemented 
the r = boundary condition by approximating the true value of ip{u,r = 0) by linear 
interpolation. Formally this is only first-order accurate. The situation is more involved, 
however: the right-hand side of @ contains an explicit factor of 1/r, which in the exact 
solution is cancelled by the very boundary condition we are trying to impose, and which 
analytically leads to g — g = 0{r) and ip — (f) = 0{r). The risk of large numerical error, 
or even numerical instability, is clear. Empirically we found that our linear approximation 
algorithm is stable, but gives rise to large errors at small u, which actually decrease with 
increasing u. We believe the reason is that at small u the scalar field is still large at r = 0, 
thus introducing a large error into the boundary condition. At large u, the field strength 
resides mostly at large r, and in this regime the standard discretization error is dominant. 
We addressed this problem, not by using higher-order interpolation, but by brute force. We 
made the grid much denser in r, typically by a factor of 16, up to and slightly beyond the 
radial scale of the initial data. This meant that the grid size was much smaller until the 
bulk of the energy in the scalar field had reached r = and then had propagated out to 
large r. In this way we achieved a reasonably small relative error everywhere, small enough 
to give us confidence in our results. 

In summary, our code is not second-order accurate, but better than first-order accurate. 
We achieved high enough accuracy to have a small relative error everywhere, but even at 
a much lower resolution the physical features we had set out to verify, QN ringing and 
power-law tails, are unambiguously present. 
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IV. SPHERICAL COLLAPSE 



In this section we present our results for the purely spherical collapse of a self-gravitating 
minimally coupled scalar field. Our expectations are based on linearized perturbations as 
analyzed in the Paper I. Because the scalar field at late times has small amplitude, it 
seems plausible that the late-time fields can be viewed, roughly, as perturbations, and that 
the analysis of Paper I applies, at least approximately. Our expectations then include the 
presence of QN oscillations at late times after a black hole has formed, and power-law tails 
at late times, both when a black hole forms and when it does not. The real and imaginary 
part of the QN frequencies - the least damped mode is the only one visible in practice - are 
fixed numbers divided by the black hole mass. The powers of the linear theory tail are fixed 
integers. When there is a static scalar field present on the initial null slice, the amplitude 
of the power-law tail, in the linearized case, is also determined. The analysis in Paper I 
suggests that these predictions for a scalar test field should hold also, approximatively, if 
the spherically symmetric scalar field evolves under the influence of its own gravitational 
pull rather than on a fixed Schwarzschild background. 

We examined two one-parameter families of initial data. In the first family the field is a 
Gaussian in r; we consider the Gaussian to represent a typical collapsing "shell of matter." 
There is a scale invariance in the problem which can be used to set the center of the Gaussian 
at r = 1; the remaining physical parameters are the Gaussian width and amplitude. A black 
hole will form from these initial data if the amplitude is sufficiently large, or if the Gaussian 
width is sufficiently small. 

As a second family of data, two different forms of (f){u = 0, r) are joined. For r less than 
some joining radius rj, the solution is constant at = 00, and for r > rj, the initial data 
is taken to be that static solution of (|l|) which is well behaved as r ^ oo (i.e., the solution 
which falls of as r~^). This solution is called "static-static" because it gives rise to a 
static spacetime in the domain of dependence of either the inner (r < rj) or outer (r > rj) 
initial data. Without loss of generality we can set the joining point at rj = 1. The only free 
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parameter is then A, the field amphtude at that point. Again this one-parameter family 
should contain spacetimes with and without an event horizon. 

Static-static data are the boundary between generic "shell" initial data and data which 
are not asjTiiptotically fiat. We found that when we replaced the exact static solution 
(given analytically as an implicit function in [y,|T|) by = const. /r) the resulting solution 
is qualitatively different: it corresponds to an infall of matter at arbitrarily late times. 

A first check of our code is the amplitude of static-static data at which a black hole first 
forms. We find 0.28 < 0o < 0.29, which agrees with the critical value 0o = 0.286 given in [|l|. 
For the Gaussian data we find 0.033 < A < 0.034, where A is the height of the Gaussian. 
These limits are stable under reduction of the grid size by a factor of 16, from 1/20 down 
to 1/320, as described in the previous section. 

A crucial feature of both Gaussian and static static data is that the amplitude (the value 
of 00 or the Gaussian height A) can be chosen either "subcritical" (no hole formation) or 
"supercritical" (a hole forms). By varying the amplitudes we found three different regimes 
for the resulting spacetimes. a) For the amplitude much greater than its critical value 
we found QN oscillations and tails qualitatively and, for the most part quantitatively, as 
predicted by linearized analysis, b) For the amplitude near its critical value we found the 
same power-law tails, but we found QN ringing to be qualitatively different in the Gaussian 
and static-static cases, c) For the amplitude smaller than its critical value (even marginally) 
we found the same power-law tails as for collapse, but no QN ringing. 

We now discuss in more detail the results for power-law tails and for QN ringing, starting 
with the former. The analysis in Paper I shows that for linear spherical perturbations tails 
should fall off in time (time for a distant observer, or Bondi time) as t]^^ in the case of generic 
data, and as for data corresponding to an initially static monopole moment. Figure 7 
shows log-log plots of at a constant values of the radius (here r = 10) as a function of 
Bondi time for different values of the amplitude of the initial Gaussian data. At late times 
clearly decays as a power of Bondi time, The measured exponents are in reasonably good 
agreement with the prediction of linearized theory: n = 2 for static static data, and n = 3 for 
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Gaussian data. These power-law tails develop, and have the same exponents, whether or not 
a hole forms, and the amplitude of the tails is a smooth function of the initial amplitude. For 
the development of a tail, hole formation is irrelevant. This independence of hole formation 
is explained by the argument in Paper I that tail development is due to backscatter at large 
r, and does not depend on the small r nature of the spacetime. 

In the case of static-static data, the linear perturbation analysis of Paper I predicts not 
only the exponent of the tail, but the amplitude. If our static static data were evolving on 
a fixed Schwarzschild background of mass M, the Paper I prediction would be that at late 
times and large radii (specifically t ^ r ^ M) the tail would have the form = —AM(j)Q/t^. 
It is unclear just how to apply this in general to the time dependent spacetimes of our 
present numerical investigations. Since the mass of the spacetime is evolving, and might 
fall to zero, what value of "M" should be used in the prediction for the amplitude of the 
tail? The picture of scattering underlying the calculations of Paper I would suggest that the 
correct M is some average, over retarded time, of the mass. The question of the appropriate 
mass is not crucial in one subset of cases: collapse models in which the initial data leads to 
the formation of a hole of mass not very different from the initial Bondi mass. In this case 
one might expect the linear fixed background prediction to apply. Our results, however, do 
not show this. In all cases, whether subcritical or supercritical, the amplitude of the tail is 
an order of magnitude less than the prediction of linearized theory (using the initial Bondi 
mass for M in the above expression). The reason for this is not yet clear. 

The results for QN ringing differ in an important way from those for tails. QN ringing 
is a late time oscillation associated with the "effective potential" governing the dynamics 
of perturbation fields outside a hole. Unlike the power-law tails, QN ringing depends on 
the small r nature of the spacetime. If a hole (or relativistically compact object) does not 
form, the generation of QN oscillations is not expected. This is verified by the results; we 
have found no QN ringing for solutions without black holes. The transition from solutions 
with QN ringing to solutions without is, however, qualitatively different for Gaussian and 
static-static data. 
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For Gaussian data, this transition is abrupt; QN ringing is clearly present in a marginally 
collapsing spacetime, and clearly absent in a marginally noncollapsing spacetime. This is 
strikingly illustrated in Fig. 8, showing the solutions evolved from Gaussian data with am- 
plitudes 0.033 and 0.034. We plot as a function of Bondi time at constant radius. Initially 
the two solutions are close, as one would expect from the closeness of their initial data. 
Later, the collapsing solution is delayed with respect to the noncollapsing one, and shows 
QN ringing. But after QN ringing has died away, and the power law tail dominates both 
solutions, they are again close. This result is compatible with the linearized analysis given 
in Paper I. Power-law tails are the result of backscatter, at large r, of the initial outgoing 
burst of waves. The two near-critical curves in Fig. 8 have almost identical outgoing bursts 
(since they correspond to very nearly the same initial data) so, at late times, the tails should 
be the same. The QN ringing, on the other hand, depends on how the solution develops at 
small r and hence is expected to be very different for the subcritical and supercritical cases. 

For static-static data, the QN ringing fades away gradually before the critical amplitude 
is reached from above. In consequence, 0(r = 10, t^) is nearly the same function everywhere 
for marginally collapsing and noncollapsing data; see Figs. 9 and 10. It may be that the 
QN ringing is masked, in the marginally collapsing case, by the non-QN fields, but the 
explanation at present is not clear. 

It is difficult to measure the QN frequencies for the spherical scalar field, as only one 
minimum and one maximum is visible. For an initial Gaussian of amplitude 0.06 for example, 
one measures a half-period of about 4.0. In this case the initial Bondi mass is 0.1699 
and the final Bondi mass still 0.1696, so that the mass is practically constant. The half- 
period predicted for spherical scalar perturbations [|1^] is 28.44M, so on the basis of linear 
perturbation theory we would expect a half period between 4.832 and 4.824, which agrees 
with what is seen in the numerical results within the rather unsatisfactory precision available. 
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V. PERTURBATIONS OF SPHERICAL COLLAPSE 



To understand better the results of the nonhnear evolution of the spherically symmetric 
scalar field, we have investigated a closely related problem. We consider a second minimally 
coupled scalar field v?, not coupled to the field, and we treat this second field perturbatively. 
That is, we ignore the contribution of (p to the stress-energy, and hence to the spacetime 
geometry. We study, then, the evolution of the perturbative ip field on the backgrounds 
generated by collapsing fields. Since these backgrounds are spherically symmetric, the 
perturbations can be decomposed into spherical harmonics which decouple. The equations 
of motion for a perturbation of fixed angular dependence p = <fY^{u, r)Y{^{9, cf)) are simply 

Di^T = Y^ig - mT - vT) + Yr^{i + ^)9^T. (12) 

pT^\[^TdT, (13) 

where g and g are still determined by the background solutions for cj), through and (^. 
The last term in (0) turns out to drive instabilities near r = for / 7^ 0. The simple 
expedient of reducing the step size Am to a value much smaller (we used 1/16) than that of 
condition (^) was found to remove the instability. (We are grateful to Jeffrey Winicour for 
suggesting this remedy.) 

It should be understood that the equations governing (p are not quite the same as those 
that would describe perturbations of the field itself. The stress-energy tensor is quadratic 
in 0, so that a perturbation 50 in would give rise to a perturbation in the stress-energy 
tensor that is first order in 50, and proportional to the "background" value of 0. It follows 
that there would be a first order perturbation in the spacetime geometry. The equations 
governing 5(f) would therefore not be the equations for a minimally coupled field on the 
background spacetime given by the background solution. (The situation is similar to 
that for electromagnetic perturbations of the Reissner-Nordstrom spacetime. Since there 
is a background electromagnetic field, the equations for electromagnetic perturbations are 
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not simply the Maxwell equations on the Reissner-Nordstrom spacetime.) Despite this, one 
should expect that there is not a great deal of difference in the late time behavior of (p 
and of 54>. Since the background 4> field becomes very small at late times, the stress-energy 
perturbations should be very small. 

The motivation for studying the (p field is that the scalar collapse models go beyond the 
models of Paper I in two distinct ways. First, the scalar collapses involve nonlinear field 
evolution, and second, they produce time dependent spacetime geometries on which the tails 
and the QN oscillations are developing. The introduction of a perturbation field allows us to 
separate these two complications. The dynamics of ip is purely linear, but takes place on the 
time dependent background spacetime. A study of the (p field has an additional attraction: 
For fixed backgrounds the exponents of the power-law tails depend on the multipole index 
/. For generic data the tails fall off as for initially static multipoles they fall off as 

^-(2i+2) gy computing different multipoles of (p we can check whether this I dependence 
applies also to time dependent backgrounds. A further technical advantage is that multipole 
modes of higher / have QN frequencies which are considerably less damped than the least 
damped I — mode. This allows us much more readily to "observe" the presence of QN 
oscillations and to measure their frequencies. 

As an example we turn again to the solution with Gaussian data of amplitude 0.06, an 
example we considered in the previous section. For the (p test field we choose a Gaussian 
of width 0.1, the same width as for the background (f) field, but we put the center of the (p 
Gaussian at 2.0 while the center of the Gaussian, which generates the background, is at 
1.0. This means that the perturbation "shell" goes in later than the background shell, giving 
the black hole time to form. For the I — cp field we again (as in section IV) find a half- 
period of about 4.0, and again find it difficult to determine the period with any accuracy, so 
that a meaningful comparison cannot be made with the fixed background prediction of 4.824 
(calculated from a final Bondi mass of 0.1696). For the I — 1 cp field, however, we could 
see 10 oscillations, and could measure a half-period of 1.811, which is within 0.7 % of the 
predicted value of 1.819. For I — 2 there are 20 oscillations, with the measured half-period 
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of 1.105, within 0.2 % of the predicted value of 1.104. QN ringing for these cases is shown 
in Fig. 11. 

We also examined I = 2 QN ringing on a background with considerable mass loss during 
the time of QN excitation. To arrange this a Gaussian "shell" of (p and a Gaussian "shell" 
of test field ip were both centered at r = 1.0. The field was chosen to have amplitude 
0.033. This gives an initial Bondi mass of 0.05822 which decreases to a final Bondi mass 
of 0.02035, with most of the mass loss taking place during 2 < < 4. For these masses, 
linearized theory predicts that the frequency of the least damped / = 2 mode should go up 
from 8.31 to 23.87. This is in excellent agreement with our numerical results, within the 
inevitable uncertainty in measuring a varying frequency. 

For power-law tails the most interesting backgrounds are those that are the most time 
dependent, the noncollapsing backgrounds. For these backgrounds results are shown in 
Fig. 12 for tails from Gaussian data with an amplitude of 0.020 and in Fig. 13 for tails from 
static-static data with an amplitude of 0.10. As these figures show, the late time behavior 
is rather accurately that of a power-law tail; the numerically determined exponents are in 
rough agreement with the predictions of linearized theory, and show the predicted increase 
of slope with multipole index /. 

VI. DISCUSSION 

In Paper I we revisited the argument for the existence of power-law tails of perturba- 
tions of Schwarzschild spacetime, in order to point out that it predicts these tails not only 
at timelike infinity, but also at null infinity and on the horizon. Reformulating the argu- 
ment yielded an important spin-off prediction: Power-law tails should form in general not 
only in perturbations of black holes. The analysis suggested that the tails should form for 
massless perturbations of any approximately spherically symmetric spacetime whose metric 
is approximately Schwarzschild, at least on an outgoing null cone of finite thickness (i.e., 
finite range of advanced time). 
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The central idea, as elaborated in Paper I, is quite simple. Radiation going out along a 
thick null cone (the "initial burst" of radiation) will be scattered back in by the spacetime 
curvature at arbitrarily large radius. The backscattered radiation then reaches a small radius 
again at arbitrarily large time, attenuated by a power of Bondi time. This backscattered 
radiation then evolves tails. The leading effect in this evolution is independent of spacetime 
curvature. The exponent of the tails is therefore the same if at late times there is, at small 
r, a star, empty space, or a black hole. 

In our numerical work described here we first set out to verify the existence of the tails, as 
well as of QN ringing, in a collapse situation. We chose the model of a spherically symmetric 
self-gravitating massless minimally coupled scalar field, because tails and QN oscillations 
seemed to be absent in the results of a previous investigation The explanation of that 
absence seems to be simply that one has to go to very late times to see the late-time 
features. We did find tails and QN ringing precisely as expected. Fortunately these features 
arise about an order of magnitude in time before our code is stopped by a diverging redshift. 

Extending the scope of previous work, we paid attention to the late-time behavior of 
solutions which are subcritical, just below or well below the margin of black hole formation, 
and we again found power-law tails, but not QN ringing. For what can be considered 
"generic" data (see Paper I), we found that QN ringing disappeared abruptly at the margin 
of collapse, while the appearance of power-law tails was continuous and unremarkable across 
the transition from subcritical to supercritical models. This last observation in itself is a very 
strong corroboration of the simple picture of tail formation set out above. Further evidence 
can be found in the dependence of the tail amplitude on the amplitude of the initial data, 
as shown in Fig. 14. The Bondi mass of the spacetime scales as the square of the amplitude 
of the initial data and we find that the tail amplitude scales as the cube. This suggests that 
the tails are scattered off the spacetime curvature only once, picking up a single factor of 
mass. 

In a second extension we considered non-spherically symmetric test fields on the dynam- 
ical backgrounds generated by collapsing scalar fields. We found that the exponent of the 
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tails varies with multipole index roughly as for fixed backgrounds. Furthermore we used the 
fact that QN ringing is less damped in the higher modes to make precise measurements of 
QN frequencies. For models in which the mass of the background was reasonably constant 
we found excellent agreement with the predictions on a fixed Schwarzschild background. For 
models with significant mass loss we found a shift in QN frequency corresponding to the 
changing mass of the spacetime. 

In summary, we found that the predictions for power-law tails of perturbations of 
Schwarzschild spacetime P] hold to reasonable approximation, even quantitatively, in a 
variety of situations to which the predictions might seem initially not to apply. 

In the interest of brevity and timeliness, the results reported here do not exhaust the 
interesting questions that might be asked. It will be particularly interesting to explore in 
more detail the behavior of power-law tails and QN ringing on the critical boundary between 
collapsing and noncoUapsing initial data. It should be said here that our code is probably 
capable only of much coarser accuracy than that of Choptuik [^, and will have to be 
modified for this purpose. On the other hand, our code was adequate for verifying, for our 
two families of initial data, one of Choptuik's crucial results: that for marginal black hole 
formation the mass of the hole depends in a universal way on the parameter of the family 
(in our case the amplitude). See Fig. 15. 

We wish to thank Jeffrey Winicour and Roberto Gomez for discussions and for making 
available their results. This work was supported in part by grant NSF PHY92-07225 and by 
research funds of the University of Utah. J. P. acknowledges hospitality and support from 
the Institute for Theoretical Physics at UCSB and the National Science Foundation grant 
PHY89-04035. 



18 



REFERENCES 
[1] R. Gomez and J. Winicour, J. Math. Phys. 33, 1445 (1992). 

[2] A. Abrahams, D. Bernstein, D. Hobill, E. Seidel, and L. Smarr, Phys. Rev. D 45, 3544 
(1992); R. F. Stark and T. Piran, Phys. Rev. Lett. 55, 891 (1985); L. Smarr in Sources 
of Gravitational Radiation, ed. L. Smarr (Cambridge University Press, Cambridge, Eng- 
land) . 

[3] R. H. Price, Phys. Rev. D 5, 2419 (1972). 

[4] E. Leaver, Phys. Rev. D 34, 384 (1986). 

[5] D. Christodoulou, Comm. Math. Phys. 105, 337 (1986). 

[6] D. S. Goldwirth and T. Piran, Phys. Rev. D 36, 3575 (1987). 

[7] E. Seidel and W-M. Suen, Phys. Rev. Lett. 69, 1845 (1992). 

[8] M. W. Choptuik, D. S. Goldwirth, and T. Piran, Class. Quant. Grav. 9, 721 (1992). 

[9] H. A. Buchdahl, Phys. Rev. 115, 1325 (1959). 
[10] E. Leaver, Ph. D. Thesis, University of Utah (1985), unpublished. 
[11] M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993). 



19 



FIGURES 



Fig. 1: The conformal diagram of the spherical collapse spacetime for a final hole mass 
Mf. Shown is our null grid in relation to the future event horizon future null infinity 
scrit+, future timelike infinity i+, and spacelike infinity ig. Null grid lines pointing to the 
top left are lines of constant v; those pointing to the top right are lines of constant u or ts- 
Null data are set on the bottom right side of the grid at u — uq. The curved lines from left 
to right are a) r = const. < 2Mf, b) r = 2M/ and c) r = const. > 2Mf. 

Fig. 2: Our null grid in coordinates t{u, r) and r. The u = const, {ts = const.) grid lines 
approach the event horizon in a finite interval of u which is an infinite interval of f^- 

Fig. 3: The convergence of (f){r — 10, ts) with decreasing grid size. On the vertical axis 
is the log of the Z2-norm of the error (compared to a very-small-grid numerical solution) ; on 
the horizontal axis is the norm of the relative grid size. Data on u = Uq for this solution is 
a Gaussian 0(r) with center at 1.0, width 0.1 and amplitude 0.06. The curves correspond 
to the following sections: a)0 < ts < 10, 6)10 < te < 20, c)20 < ts < SO d)30 < te < 
40 e)40 <tB < 50. The slope of a is 1.34 ± 0.07 and of d is 1.72 ± 0.04. 

Fig. 4: The convergence of 4>{r,tB = 50) with decreasing grid size. The axes and the 
initial data are the same as for Fig. 3. From top to bottom are < r < 10, 10 < r < 20, 
20 < r < 30, 30 < r < 40 and 40 < r < 50. The slope of the top graph is 1.21 ± 0.04 and 
the slope of the bottom graph is 1.70 ± 0.04. 

Fig. 5: Part of a plot of 0(r = 10) versus t^, in a region where is dominated by QN 
ringing. From bottom to top are shown runs with initial radial grid sizes 1/20, 1/40, 1/80, 
1/160 and 1/320. Initial data is the same as in Fig. 3, which corresponds to the formation 
of a black hole with negligible radiation of mass. 
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Fig. 6: Results of the same runs as in Fig. 5, but here in a region where is dominated by 
the power-law tail. Note that the higher the precision, the smaller the value of ts where the 
run stops due to an overflow, here at 53 and 57. From bottom to top the curves correspond 
to initial radial grid sizes of 1/20, 1/320, 1/160, 1/80, 1/40. 

Fig. 7: Log- log plots of 0(r = 10, is) for Gaussian data with different amplitudes. Each 
plot starts after the last change of sign of 0. The amphtudes of the initial Gaussian are 
a) 0.06, b) 0.034 (marginally collapsing), c) 0.033 (marginally noncoUapsing) , d) 0.01. The 
power law nature of all curves is clearly visible; the exponents are a) —2.74, b) —2.63, c) 
—2.63, d) —2.68, compared to a hnearized theory prediction of —3. Only the two collapsing 
cases, show QN ringing. 

Fig. 8: A closeup on 0(r = 10, t^) in the region of QN ringing for initial Gaussian data. 
Solid hne: Amplitude 0.034, which collapses marginally, showing QN ringing (the small 
feature around t ^ 5). Dotted line: Amplitude 0.033, which marginally does not collapse. 
Note the complete absence of QN ringing. 

Fig. 9: A closeup of 0(r = 10, t^) in the region of QN ringing, for static-static data a) 
Amplitude 0.50, which collapses, and here shows QN ringing, b) Amplitude 0.35, showing 
only faint QN ringing, c) Amphtude 0.29, which marginally collapses, d) Amphtude 0.28 
which marginally does not collapse. There is no qualitative difference between c) and d), in 
contrast to the case of Gaussian data depicted in Fig. 8. 
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Fig. 10: Log-log plots of 0(r = 10, ^b) for static-static data with different initial ampli- 
tudes 00- Each plot starts after the last change of sign of 0. The amplitudes are a) 0.5, 
6) 0.29 (marginally collapsing), c) 0.28 (marginally noncoUapsing) , d) 0.01. The power law 
nature of all curves is clearly visible; the exponents are a) —2.08, b) —1.98, c) —1.95, d) 
—1.86, compared to a linearized theory prediction of —2. The power laws are a) —2.08, b) 
—1.98, c) —1.95, d) —1.86, compared to a linearized theory prediction of —2. Only curve a 
shows clear QN ringing. 

Fig. 11: Test fields </?P(r = 10, ts) on a Gaussian background, amplitude 0.06. The three 
graphs are I — 0, I — 1, I — 2 {in order of increasing frequency). 

Fig. 12: Log-log plots of test fields <^5"(r ~ 10, is) for different multipole indices / on 
one background spacetime. The background is evolved from Gaussian initial data for 4> 
with amphtude 0.02 (noncoUapsing). The initial data for the test field is in each case a 
Gaussian (the test field amplitude, 1.0, has no significance). From top to bottom the curves 
correspond to / = 0,1,2,3. The best fit for the power law exponents are —2.77, —3.95, 
—5.94, —8.34, compared to predictions of —3, —5, —7 and —9. 

Fig. 13: Log-log plots of test fields (pi^{r — 10, ts) for different multipole indices I on 
one background spacetime. The background is evolved from static-static data for (p with 
amplitude 0.1 (noncoUapsing). The initial data for the test field is in each case a Gaussian 
(the test field amplitude, 1.0, has no significance). From top to bottom I — 0, then I — 1, 
I — 2 and I — 3. The best fit for the power laws exponents are —2.70, —3.66, —5.52, —7.26, 
compared to predictions of —3, —5, —7 and —9. 
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Fig. 14: Log-log plots of the initial Bondi mass and the amplitude of the power-law 
tail (j){r = 10, ts), both versus the amplitude of initial data. The upper graph of each pair 
(full dots) represents Gaussian data, the other (empty dots) static-static data The region of 
dense data points marks the collapse-noncollapse transition on each graph. Both the initial 
Bondi mass and the amplitude of the tail scale rather precisely as powers of the initial data 
amplitude for noncollapsing data. The powers are, for the mass 1.99 (Gaussian) and 1.97 
(static-static), and for the tail 3.22 (Gaussian) and 3.01 (static-static). 

Fig. 15: Mass of black hole formed versus difference between the amplitude of the initial 
data and the critical amplitude. On the vertical axis, In(mass). On the horizontal axis 
ln[{p — p*)/p*], where p is 0o or A and its critical value. Empty squares denote static- 
static data, full squares denote Gaussian data. The fact that the static-static data points are 
further to the left means that they were measured closer to criticality. This was necessary 
because the universal power-law behaviour of the mass develops only closer to criticality 
for static-static than for Gaussian data. The deviation from this behaviour is clear in the 
plot for the empty squares for (0o — 0o*)/0o* > For smaller values the empty squares 
fall approximately on a straight line. The sets were shifted vertically (but not horizontally 
as was done in Ref. Jill) with respect to each other in order to place them on one line. 
This corresponds to an overall constant multiplicative factor in the black hole mass which 
is not universal but depends upon the family of initial data (Gaussian or static-static). We 
assumed = 0.03280 and 0o* = 0.2860. Formally, the best fit to the slope (power-law) 
for Gaussian data alone is 0.39, for static-static data alone 0.31, and combining both sets 
0.35. A full analysis of the uncertainties in these slopes has not been carried out, but we 
assume that within numerical accuracy the slopes agree with each other and with the slope 
0.37 found by Choptuik [jlT . 
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